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RELATIONS BETWEEN TRANSFER MATRICES AND NUMERICAL 
STABILITY ANALYSIS TO AVOID THE nO PROBLEM * 

R. PEREZ-ALVAREZ t, R. PERNAS-SALOMON AND V. R. VELASCO § 

Abstract. The transfer matrix method is usually employed to study problems described by 
N equations of matrix Sturm-Liouville (MSL) kind. In some cases a numerical degradation (the 
so called Qd problem) appears thus impairing the performance of the method. We present here a 
procedure that can overcome this problem in the case of multilayer systems having piecewise constant 
coefficients. This is performed by studying the relations between the associated transfer matrix (T) 
and other transfer matrix variants. In this way it was possible to obtain the matrices which can 
overcome the Qd problem in the general case and then in problems which are particular cases of 
the general one. In this framework different strategies are put forward to solve different boundary 
condition problems by means of these numerically stable matrices. Numerical and analytic examples 
are presented to show that these stable variants are more adequate than other matrix methods to 
overcome the Qd problem. Due to the ubiquity of the MSL system, these results can be applied to 
the study of many elementary excitations in multilayer structures. 

Key words. Transfer matrix, matrix Sturm-Liouville problem, numerical stability, quadratic 
eigenvalues, Qd problem 
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1. Introduction. The study of elementary excitations in multilayer systems 
(heterostructures) continues to be a very active field of research due to the multiple 
applications of these systems for the design of devices with composite materials. In re¬ 
cent years magneto-electro-elastic materials [1] and piezoelectric multilayer structures 
[2, 3], among other systems, have been the object of many studies. The associated 
transfer matrix method [4] T is one of the theoretical techniques most employed in 
the study of these systems. From the formal point of view this method is very ade¬ 
quate for the study of problems related with multilayer systems. It reflects in a very 
simple way the linearity of the problem, based on the fact that any solution can be 
expressed as a linear combination of a chosen basis of the corresponding functional 
space [4]. On the other hand, for several practical applications to different problems, 
the method is hampered by numerical instabilities, the most common one being called 
the rid problem [4]. The name associated to this numerical instability derives from 
the elastic waves studies where this instability is present at high frequencies n and/or 
big thicknesses (d) of the layers. 

A description of this problem was given in [5] when studying wave propagation in 
layered elastic media at high frequencies. In this work the origin of the problem was 
assigned to the large frequency-thickness (/d) products. It was found that the modal 
calculations presented numerical difficulties and a matrix formulation, the d-matrix 
method, was proposed to deal with them. Another approach closely related with the 
scattering matrix method is the reflectivity matrix method [6]. Nevertheless the high 
frequency-thickness product instability has been a persistent feature in the study of 
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wave propagation in layered media as it can be seen in a later review [7]. The name 
fid problem was coined in [4]. There it was noticed that the expressions producing 
this instability are of the form fid where fl stands for a function of the frequency or 
the equivalent magnitude for other elementary excitations. 

Different techniques have been developed to deal with this problem. Some of 
them, as the global transfer matrix [4, 3], involve matrices with dimensions increasing 
with the number of layers forming the system. It is then clear that for systems 
including many layers the method will require big amounts of computer memory and 
time. 

Other approaches employ transfer matrices with dimensions independent of the 
number of layers. Among them we can find the Stiffness matrix method [8, 9, 10, 11] 
(E), the Scattering matrix method [11, 12, 13] (S) and the method of the hybrid 
compliance-stiffness (or simply hybrid matrix) [11, 14, 15] (H). These methods have 
been mainly used in studies of elastic waves propagation in anisotropic systems and 
acoustic waves in piezoelectric systems. 

All these studies have been performed in a separate way. There is no clear picture 
of the usefulness and limitations of the different approaches. Our aim is to give an 
unified view of the problem and present the most adequate transfer matrix variant 
to solve different problems. To this end we shall consider a general system of N 
differential equations of the matrix Sturm-Liouville (MSL) kind [4]. In this general 
framework the study of the expressions relating each different matrix with T will 
allow to understand how the different matrices elude the T numerical instabilities. 
To extend the use of these transfer matrix variants H, E and S' to a wider range of 
physical problems involving multilayer systems we shall study the numerical stability 
of each matrix variant. In addition we shall present different strategies to be used in 
the case of common boundary value problems as superlattices or finite sandwiches in 
terms of H, E or S. 

We must stress that our procedure can overcome this problem in the case of mul¬ 
tilayer systems having piecewise constant coefficients. This is an important problem 
and covers many cases of practical interest. 

Among the big amount of work done on this problem using other methods we 
can mention those based on the sextic formalism for the linear elasticity [16]. In this 
scheme the matricant matrix was introduced together with the impedance matrix 
and the two point impedance matrix [17]. This approach allows to deal with systems 
having inhomogeneous coefficients. Stable methods to compute the matricant and the 
impedance matrix with special integration schemes [18] and an alternative method 
based on the resolvent of a propagator have been presented recently [19, 20]. In 
these works the chain rule for the resolvent, together with a differential equation of 
Riccati kind for the obtention of the resolvent in continuous inhomogeneous media 
are presented. The resolvent is well adapted to get the spectrum and fields in these 
systems. 

The general character of our approach allows the extension of the transfer matrix 
variants use to problems whose systems of equations are particular cases of the MSL. 
We shall illustrate, for example, the hybrid matrix numerical stability by numerical 
studies of the shear horizontal surface waves in piezoelectric multilayer systems. 

In Section 2 we present the master equation of the matrix Sturm-Liouville system 
of equations. In Section 2.1 we introduce the quadratic eigenvalues problem leading to 
the linearly independent solutions of the system together with their eigenvalues for a 
homogeneous medium. We define in Section 2.2 the associated transfer matrix T and 


TRANSFER MATRIX NUMERICAL STABILITY AND THE fid PROBLEM 


3 


introduce the form employed in the analysis of the numerical stability of the variants 
H, E and S. Section 2.3 introduces the fid problem together with the analysis of the 
T characteristics that can be the source of this numerical instability. Afterwards, the 
numerical stability of the hybrid matrix, Section 3.1, and the stiffness matrix. Section 
3.2, are studied through their respective relations with T in an homogeneous domain. 
The analysis for the scattering matrix is presented in (Section 4.1). The composition 
rules for the different matrices considered here are analyzed in Section 4.2. Section 5 
presents the strategies to solve several boundary problems in terms of H, E or S. A 
numerical example demonstrating the numerical stability of the hybrid matrix is also 
presented together with an analytic study of the well known Kronig-Penney model. 
Conclusions are presented in Section 6. 

2. Matrix Sturm-Liouville system of equations (MLS). A matrix Sturm- 
Liouville problem emerges naturally in a wide range of physical and technological 
problems (see, for example Refs. [4, 21, 22], and citations therein). In this wide 
range of problems there are many belonging to the elasticity theory (see for example 
[23]), electromagnetism [24] and several other areas of classical physics. Some of these 
problems can be quite complicated as the magneto-electro-elastic waves [25]. A matrix 
Sturm-Liouville problem appears also in Quantum Mechanics and Solid State Physics. 
Particularly the Envelope Function Approximation (EFA) [26, 27] generates a massive 
class of systems of equations that follow the Sturm-Liouville equation in matrix form. 
Initially many of these systems of equations are three-dimensional, but in layered 
systems, as outlined in Figure 1, the normal modes can be chosen as exponential of 
in ■ p multiplied by some function of the variable 2 , the coordinate perpendicular to 
the interfaces. We denote by p = xCx + ycy the position vector in the plane of the 
interfaces and by k = k^Cx -I- KyCy the corresponding wavevector. In this way the 
equations of motion take the Sturm-Liouville form, namely: 


( 2 . 1 ) 


A. 

dz 


B{z) 


dF{z) 

dz 


+ P{z)-F{z) 


+ Y{z) 


dF{z) 

dz 


W{z)-F{z) = 0 


This defines the matrix differential operator L(z).The unknown F( 2 ;) is the field 
under study: electronic wavefunctions, or envelope functions, if we deal with elec¬ 
tronic states, vibration amplitude for elastic waves, or components of the electric field 
in some electrodynamic situations. In the case of the Full Phenomenological Model 
(FPM) for polar optical modes in heterostructures [21] the unknown field has several 
components: three mechanical amplitudes and a component which is interpreted as 
a coupled electrostatic potential [21]. The coefficients B{z), P{z), Y{z), and W{z) 
are square matrices of order being N the number of coupled second order differ¬ 
ential equations forming the system (2.1). These coefficients characterize the physical 
properties of the materials forming the multilayer system: dielectric constants, elastic 
coefficients, etc. As the multilayer structures studied here involve different materials 
these coefficients will be different for the different materials. The dot • means standard 
matrix product. 

As the linear differential form is defined from (2.1) 


(2.2) a{z) = B{z)-^ + P{z)-F{z), 

then the first integration from z — e to z + e shows that A( 2 ;) is continuous for every z 
along the multilayer structure. The continuity of F(z) and A(z) along the structure 
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allows to obtain the composition rule for the transfer matrices defined from these 
magnitudes. 

2.1. LI solutions. Quadratic Eigenvalues Problem. In the case of an ho¬ 
mogeneous medium the differential equations system ( 2 . 1 ) takes the following form 

(2.3) B ■ F"{z) + {P + Y)- F'{z) + W ■ F{z) = 0 . 

In this simple case the linearly independent (LI) solutions of the differential equations 
system (2.1) can be expressed by means of exponentials [28, 29] 


(2.4) 


F(^) = . 


The eigenvalues k are obtained from the zeros of the secular matrix determinant: 


(2.5) &{k) = -k^B + ik{P + Y) + W . 

Now we are dealing with a quadratic eigenvalues problem (QEP) [22]. If matrix 
B is regular (Det[J5] 7 ^ 0) we have a set of eigenvalues K = {kj,j = 1,2, ••• ,2N} 
and the corresponding eigenfunctions Fj{z) = F^o exp[ikjz]. The amplitudes Fjo 
multiplied by a constant are obtained from the homogeneous linear equations system: 


(2.6) ©(fc,) • F,o = 0 . 

The multiplicative constant is usually obtained by a normalization condition. 

In the following, we shall always assume B^ = B, = W and Y = —P\ 
in order to ensure formal hermiticity of the operator L(z), see Ref. [4]. In this case 
the eigenvalues of the QEP satisfy the general property of being real or appearing in 
pairs: kj and its complex conjugate k*. 

2.2. Associated Transfer Matrix for the MSL equations system. We 

shall define the associated transfer matrix T{a : z, ^o), which transfers the amplitudes 
F and the linear differential form A in a domain a, as in [4]: 


(2.7) 


F(a : z) 
A(a : z) 


T{a : z,zq) 


F(a : zq) 
A{a : zo) 


From now on we shall suppress the zonal argument a. Following the algebraic 
and analytic methods to calculate the matrix T(z,zo), given in [4], we shall have: 


( 2 . 8 ) T{z,zo) =Q{z)-Q{zo)-\ 

where the auxiliary matrix Q{z) is formed by a basis of eigenfunctions Fj(z) and of 
the linear differential forms Aj{z) = B 

dz 


Q{z) = 


(2.9) 


Fi(z) F2(2:) ... F27 v(-z) 

Ai( 2:) A2(2:) ... A2Ar(-z) 
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For an homogeneous domain a, with constant B, P and W, we can choose the 
eigenfunctions Fj{z) = Fjo and after some manipulations on (2.9), we can 

separate the factors F^o from the exponentials in the form: 


FOiv 

^02Ar 


O 

Kf 

1 

G 


(-^o) 0 

■^On 

-^02Ar 


0 Flk^^{z - zo) 


0 Uk^^izo) 


The submatrices Flfcjy (d) and {d) are diagonal an the jth element is the 
exponential Ao^, Fo2 „ and are square matrices of order N whose 

elements are obtained in terms of the constant F^o and the corresponding A^o- In 
our notation the subindex {N} denotes that j = 1, 2, • • • ,N and the subindex {2A} 
means that j = A + 1, A + 2, • • • , 2 A. 

By substituting Q{z) in (2.8) and considering d=z — zq we have: 


(2.10)T(d) 


^Oat ^02iV 


o 

G 



■^02N 


0 nfe2iv (d) 


Aqat ■^02N 


The matrix T(d) appearing in (2.10) can be interpreted as the associated transfer 
matrix (ATM) relating the vector [F(z) A{z)^ in the boundaries of an homogeneous 
domain with thickness d. 

As the linear form A.{z) is continuous along the interface separating two adjacent 
domains, the ATM has the chain property. Then for an ensemble of /j, layers sketched 
in Figure 1, the system ATM is obtained from the following matrix product: 


(2.11) T{Zr, Zi) = T{Zr - Zf^-i) ... T{z2 - Zi) ■ T{zi - Ze), 

where Zi, zi, Z 2 t ■ ■, Zr are coordinates of the interfaces matching the different domains 
of the multilayer structure. 

We shall start now the study of the T characteristics which can be the source of 
the fid problem in the numerical calculations. With this knowledge we shall study 
later the numerical stability of the iT, E and S matrices, by means of their relations 
with T. 

2.3. fid problem. From (2.10) we can obtain expressions for the analysis of 
the numerical instability of the T matrix elements for any A. For real eigenvalues 
(allowed regions) we have: 


2N 

(2.12) T,, = 2 Aisj [cos(A:j d) + i sin(fcj d )], 

1=1 

whereas for complex eigenvalues (forbidden regions) we shall have combinations of 
decreasing and increasing exponentials: 

N 

(2.13) Ti, = ^ Cisj T, ^ (l ± . 

1=1 


The coefficients Aisj, Cisj and Disj are expressed in terms of the elements of Fq^^, 
Aq^j,, Foaw and Aqj^,. In (2.13) we have separated the kj eigenvalue real 'Si{kj) and 
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m 


m + 1 



R 


Zi Zi Z2 Z3 ■■■ Zm-1 Zm Zm + l ' ’ ’ 2^-3 Z^-2 Z^-i Zr Z 


Fig. 1. General scheme of the system under study. The system consists of layers sandwiched 
by two semi-infinite external domains: L (left) and R (right). According to our convention, the 
layer m is bounded between interfaces (m — 1) and m with coordinates Zm—i o^^d Zm respectively. 

imaginary parts. The real part is included in the Tj = factor having a 

bounded value. 

For real eigenvalues the T elements are represented by means of trigonometric 
functions which are bounded by ±1. In this case the fid problem does not appear 
when the product (kj d) increases. On the other hand, for complex eigenvalues the 
mixing of terms with increasing and decreasing exponential values present in (2.13) 
may give rise to this nuirrerical instability. 

For increasing kj d leading to ^ ^ u {unit roundoff) the (l ± ‘^) 

operation is rounded to 1.0 by the computer. Thus the result (1.0) will have a round¬ 
off error. The number u {unit roundoff) is the machine precision, that is, the value 
to be added to 1.0 to produce a result different from 1.0. This number can be cal¬ 
culated as It = 5/3^“* [30], where /3 is the base of the floating point number system 
and t its precision (can be understood as the number of digits used to give a value). 

The roundoff error is defined as the difference between the calculated approximation 
of a number and its exact mathematical value. When the roundoff result is 1.0 the 
absolute value of this error Eats is bounded Eabs ^ u [30]. In the double precision 
decimal system (/? = 10, t = 16) we have u = 5 x 10“^®. 

If we assume that the calculation of a term Cisj tj (l ± '^) 

of (2.13) is performed with roundoff, then the result will be affected by a roundoff 
error with absolute value Er given by: 


(2.14) 


Ej‘ ^ C^isj 'Zj 




Depending on the numerical problem under study the right-hand side in (2.14) can 
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have a high value and also a big Ej. error. The roundoff error can be accumulated when 
the final result to be obtained (e.g., eigenvalues or parameters of a given problem) is 
preceded by a sequence of calculations prone to roundoff errors. In these cases the 
error can dominate the calculations thus giving a very inaccurate final result. When 
this happens we are in the presence of the numerical instability called Od problem. 

In practice it is quite easy to deal with problems in which the T determinant is 
constant and equal to one. This can be used as a test of the numerical accuracy in the 
real calculations. When the numerical instability is present the Det[T] takes values 
quite different from the exact one, being in some cases several orders of magnitude 
bigger or smaller than 1.0. 

The expression (2.13) shows also clearly that the Tig elements increase indefi¬ 
nitely when exponential argument |3(A:j)| d —> oo. In this case the T matrix overflows 
and cannot be calculated numerically. Thus it is clear that in the T numerical appli¬ 
cations we can find two kinds of numerical instability: the fid problem and the matrix 
overflow. 

3. Hybrid matrix and Stiffness matrix of the MSL system. We can define 
new matrices in the domain a, where T was defined, by changing the arrangement 
of the F(2:), A(2 :), F(zo) and A{zo) vectors in (2.7). Some examples are the Hybrid 
Compliance-Stiffness matrix (H) and the Stiffness matrix (E): 


(3.1) 


(3.2) 


F(a : 
A(q; 

Zq) 

:z) 

II 

z;zo) ■ 

A(a 

F(a 

Zo) 

:z) 

A(q: 

A(q: 

zo) 

:z) 

= E{a 

z;zo) ■ 

F(a 

F(a 

Zo) 

■■z) 


The Hybrid Compliance-Stiffnes matrix was employed in Ref. [14] as a stable 
variant to study the propagation of an acoustic wave in an anisotropic multilayer 
system. The acoustic wave equations of motion are a particular case of the system 

(2.1) including the displacement vector as F(2) and the force vector normal to the 
interfaces as A(z). 

Following this procedure we can define up to 24 interrelated matrices related 
among them. In fact we obtain 12 different matrices and their respective inverses. 
Among them we find H~^ and E~^. The matrix £^~^is known as Compliance 
matrix, see Refs. [14, 8]. By taking as reference the expressions defining T, H, E, 
H-^ and E-^ is possible to obtain from them other three different matrices 
which will exhibit a similar numerical behaviour. A first matrix is obtained by per¬ 
muting among them the positions of the vectors in the right-hand side of the matrix 
taken as reference (e.g., the A(a : zq) and F(q; : z) vectors in the right-hand side of 

(3.1) . The second matrix is obtained by means of the former operation applied to the 
vectors in the left-hand side of the matrix taken as reference. The third one is the re¬ 
sult of both permutations. The Appendix A shows, by means of the relations between 
the matrices, that the matrices defined in this way will have a similar behaviour from 
the numerical point of view. 

We obtain T{—d) by inversion of (2.10). Thus will have the same numerical 
behaviour than T. When inverting the expressions for H{d) and E{d) the result is 
the permutation of the Fsubmatrix with the Ao;,^ one and of the Fsubmatrix 
with the Ao 2„ one. Thus and E~^ will have also a numerical behaviour similar 
to those of their counterparts. 
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3.1. Analysis of the numerical instability of the Hybrid matrix of the 
MSL system. The following relations can be obtained from Eqs. (2.7) and (3.1): 


(3.3) 


H = 


-[Tn]-i.ri2 

r22-T2i-[Tn]^i.Ti 


[Ti 


1-1 


On the other hand, equations (2.7) and (2.8) exhibit an important property. 
Equation (2.8) leads to a unique ATM independently of the chosen LI solutions base. 
As a consequence the hybrid matrix obtained from the relations (3.3) will be in¬ 
dependent also from the solutions base chosen to build up T. Then, for simplic¬ 
ity, we consider that the ATM expression (2.10) was built from a base of solutions 
Fj(z) = in such a way that nfc;v(d) contains the eigenvalues with posi¬ 

tive imaginary part 3(fci, k 2 , ■ ■ ■ kpf) > 0 and nfc 2 j^((i) the eigenvalues with negative 
imaginary part S(/cAr+i, kiq+ 2 , ■ ■ ■ ^ 2 ^) < 0. In this way the submatrices H^j, (d) and 
nfc 2 „(—d) reduce to the order N nil matrix {On) when the thickness d —> oo whereas 
the elements of nfc„(—d) and nfe 2 j^(d) tend to infinity. 

Appendix B contains the expressions for the N order partitions: Tn, T 12 , T 21 
and T 22 obtained from (2.10).With the help of (3.3) we have: 


(3.4) 


Hii — [A 02 JY • Fq^^ 721 • Hfe^( d) ■ Fg^ • Fqj^ ■ Tlk2N{d) ■ 7i2 ] 


[A, 


On 


On 


■ 722 • nfc 2 ^(-d) • F, 


1 Fo, •nfc„(d)-7ri'] ' 


02iV 


(3.5) 


ri —1 


Hi2 - 712 • nfc2„(-d) • Fp^^ 

• [In - Fo„ • nfc„(d) ■ Ag^ • Aq^^ • Uk^Ni-d) ■ 


(3.6) 


H21 — [Aq;^ — H22 • Fo„] ■ nfc„(d) • ^21 + 
[■^02iv H 22 ■ Fqjjv] ■ Flfcjjv {d) ■ 722 


[Fo„- 

■^On F 02 N 

FIfc 2 N {d) 

■ A~^ 

^ 02 N 

* ^On 

■ Flfc^ (- 

-d). 

Aoi]"V 

[F 02 JV 

' ■^ 02 N' ~ 

• Flfc;,, (d) 

■ A-^ • 

^02N 

■ FIfe 2 JV ( 

-d) 



The coefficients 711 , 712 , 721 and 722 are obtained in terms of Fo;„, Fqj^, Aq;,, 
and Ao 2 „ as indicated in Appendix B. 

When d increases indefinitely the expressions (3.4-3. 6 ) are reduced to: 


(3.8) 

Hn 

d—>-ao 

= Fo„ ■ Ag^ • [In - Oat] ^ = 

Fojv • A 

-1 

On 

(3.9) 

Hi2 

d—>-cc 

= Oat • [Jjv — Oat] ^ = Oat 




H 21 

d—>-(X) 


^02N -^02N ' ^ 

1 

0 

1 

• Fqjv 

(3.10) 



+ 

^02N ^ -^02^7 [-^A^ 

•nfc2„(d)| 

(3.11) 

H 22 

d—>-(X) 

= Ao2„ ■ Fp^’^ • [In - Oat] ^ 

-^027^ 

■p—1 

^ 02 N 
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We denote by [/jv — Ojv] the N order identity matrix obtained with roundoff, 
whose elements are characterized by a roundoff error with absolute value Er ^ u 
{unit roundojf). Thus, the H matrix elements will be characterized by an error whose 
absolute value is of order of u. These results show that the MLS matrix H converges 
to finite values, without significant precision loss, when d increases indefinitely. 

On the other hand when d 0 we obtain immediately from (2.10) that T = I 2 N 
and by substituting its partitions of order N in (3.3) we obtain: 


(3.12) 


= 


Oiv In 
In Ojv 


thus H also converges in a numerically stable way when d —> 0. 

3.2. Numerical stability of the Stiffness matrix of the MSL system. 

From the expressions (2.7) and (3.2) we derive the following relations: 


(3.13) 


T21 - T22 • [Ti2]-i • Tn T22 ■ [Ti2]-i 


The Stifness matrix obtained from equation (3.13) will be, as the H matrix, 
independent of the base of the LI solutions chosen to build T. Because of this we 
consider also the ATM coming from the expression (2.10), which was obtained from 
a base of solutions Fj ( 2 :) = where nfcj^(d) contains the eigenvalues with 

positive imaginary part: 3(fci, ^ 2 ,... ^at) > 0 and nfc 2 jy(d) contains the eigenvalues 
with negative imaginary part: '^{kN+i, kN+ 2 , ■ ■ ■ k 2 N) < 0. 

Following the same procedure employed for H we substitute in (3.13) the expres¬ 
sion of the partitions Tn, T 12 , T 21 and T 22 given in the Appendix B and calculate 
the limit of the partitions of E when d —> 00, to obtain: 


(3.14) 

F'li 

d—>CX) 

= • [In - Oat] ^ = 

Aon 

•Fo'^ 

(3.15) 

Ei 2 

d—>CX) 

0 

II 

7 

0 

1 

0 

II 




E21 

d—>-cX) 


‘-Oat -^02Ar * ^02iv ’ ^ 

- Oat] 

• Fo.] 

(3.16) 



+ 

■^02Ar ■^02iv [-^A^ 

-]■ 

FlfcjN- (rf)l 

(3.17) 

E22 

d—>CX) 

= Ao 2„ • Fg^^^ • [In - Oat] ^ 

-^02N- ■ Fg^jy 


These results show that the E matrix also converges to finite values without a 
significant precision loss when d grows indefinitely. On the other hand, when d —>■ 0, 
we know that T = I 2 N, which means that T 12 = T 21 = On and then the E matrix is 
not numerically computable (overflow) as is directly obtained from the relations (3.13). 
Let us now assume that d is very small but not enough to provoke the overflow state. 
From (3.13) we can express the partition E 21 in the following form: 


(3.18) E 21 = -(In- T 21 ■ T^i^ • Ti 2 • T 22 ^) • T 22 • T)- 2 ^ • Tn, 

then for a sufficiently small d this partition will be the object of the roundoff in the 
first place, giving: 
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(3.19) E2i\d^o = — [In — Oat] • T’i2^. 

Unlike the limits given in (3.14)-(3.16) the term [In — Oat] subjected to the round¬ 
off, multiplies now a term T ^2 whose value can be big enough to affect the Stiffness 
matrix due to the roundoff error and then gives rise to the fid problem. 

4. Scattering Matrix and Coefficients Transfer Matrix for the MSL 
system. In the case of the Scattering Matrix (S) its relation with T is not a direct 
one (because there are other matrices involved) as in the relations studied previously. 
A possible way to relate S with T is by using the Coefficients Transfer Matrix (K). 
We need to use the direct relation S-K and the indirect one K-T. Being known a 
base of solutions Fj(Q;, z) in a domain a the general solution F(q;, z) of the differential 
system (2.1) can be written as: 


2N 

(4.1) F{a,z) =Yjaj{a)Fj{a,z) . 

3 

Let be a+/“ (a) the iV-vector formed by the coefficients aj{a) of the amplitudes 
travelling to the right/left. Then we shall denote as If (R, L) the Coefficients Transfer 
Matrix transferring the ensemble of coefficients sl^/~ from domain L to domain R: 


(4.2) 


a+(R) 

a-(R) 


K{R,L) 


a+(L) 
a (L) 


The term Scattering Matrix is widely used in the literature and can be defined in 
different ways. Here we shall use the definition and notation iS'(R; L) employed in [4]: 


(4.3) 


a-(L) 

a+(R) 


S'(R;L) 


a+(L) 

a-(R) 


From these definitions we obtain a direct relation between S and K: 


«=r -[if22]-'-if21 [if 22]-' 

^ ’ [Ki,- If 12 • [if 22 ]-' • if 21 if 12 • [if 22 ]-' • 

By taking into account that between the domains R and L there is an intermediate 
region M (can be a single or a multiple layer) described by a T matrix, it is possible 
to obtain [4]: 


(4.5) if (R, L) = [Q(R : z.)]"' ■ T(z„ zi) ■ Q{L : zi) , 

where Zt/i are the coordinates of the interfaces matching the intermediate region M 
with the external domains R (to the right)/L (to the left). The matrix Q{z) for an 
arbitrary domain is given in (2.9). Expression (4.5) shows clearly that the matrix K 
and consequently S depends on the base of LI solutions chosen to build the matrix 
Q. It is a common practice to choose a reduced base in z^/i^ that is a base tending 
to unity in z^jg. 
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4.1. Analysis of the numerical stability of the Scattering Matrix (SM). 

In the first place we substitute in (4.5) the expressions (B.2) giving the T partitions 
when d —> 00 . Now we substitute in (4.4) the expressions obtained for the partitions 
of IC and obtain: 


5^11 U^oo = — [7 i2^ ■ (-^N + Otv) • Q(z£)i2 + 722^ ' JV + ^n) ’ Q(z£)22] 
(4.6) ■ \pii2 ■ {In + Oat) • Q{zi)ii + 722 ^ • {In + O^v) • Q{z£)2i] ; 



‘S'12 d—»-GO “ 

II 

TT 

■f Oat) • Q{zi)i2 + 722^ ■ 

{In + On) ■ 

Q{z£)22] 



'nfc2N'( d,)\ 

d—>00 ' )21 * Fo2;v + {Qi^r) 

~^)22 ■ Ao2„] 

(4.7) 


= Oat ; 





^21 d—»-GO “ 

= [(Q(^r)-') 

11 ■ F 02 N + {Q{Zr) ^)l2 

' Ao 2N'] ' Flfcj^ (d) d—VCX) 



• [ 71 - 2 ' • {In 

+ Oat) • Q{z£)ii + ')22 

* (In + Oat) 

• Q{z£)2i\ 

(4.8) 


—Identical 

= Oat ; 




5^22 d-»GO = 

= [(Q(^r)-') 

11 ■ F 02 N + {Q{Zr) ^)l2 

* -^ 02 ^ 7 ] 


(4.9) 


■ [{Q{Zr)~^ 

)21 ■ F 02 N + {Q{Zr) ^) 22 -Ao 2 ;v] 



We have used the notation Q{zr) instead of Q(R : Zr) and Q{ze) instead of 
Q{Ij : z^) to simplify the expressions. The term [Jat + Oat] is an identity matrix of 
order N obtained by roundoff whose elements have a roundoff error with absolute 
value Er ^ u {unit roundoff). Because of this the matrix elements of S will have an 
error with absolute value of the order of u. These results show that the S matrix of 
a MSB system converges also to finite values without significant precision loss when 
d increases indefinitely. 

On the other hand, if the intermediate region M thickness is nil (no M region, 
zi = Zr = Zs) we obtain from (4.5) that X(R, L) = [Q(R : -2s)]“^ ■ Q{L : Zg) and we 
can obtain S without trouble. This means that the S matrix of the MSL can avoid 
the Od problem and converge in a stable numerical way when d ^ 0 . 

4.2. Composition rules. The hybrid matrix relating the field and the 

linear form in the positions Zm-i and Zr of the sketch shown in Figure 1 can be 
described as the hybrid matrix of the structure formed by the layers m, m + 1 ,... ,/i: 


(4.10) 


F(m : Zm-i) 

— Ff(™) . 

> 

1 

A{n : Zr) 

— JfJ. 

F(/r : Zr) 


We use here the supraindex m among parentheses to denote the hybrid matrix 
of the structure being considered in a similar way to that employed in Ref. [14]. 
The partitions of the matrix given in (4.10) can be expressed in terms of the 
matrix partitions corresponding to the structure including the layers from m + 1 to 
fj, and of the matrix iT"* given by (3.1) relating the field and the linear form in the 
layer m borders zq = Zm-i and z = Zm- We must take into account the continuity 
of the field and the associated linear form in Zm, F/A(to + 1 : Zm) = F/A(m : Zm)- 
Then we obtain the following composition rule: 






12 


R. PEREZ-ALVAREZ, R. PERNAS-SALOMON AND V. R. VELASCO 


Trim) _ Trm , tt 

■tl — ±1 ~r -tl 


12 ■ -f^ll 


(m+1) 


■jnr{m) _ rrm 
^12 — ^12 


In^H 


(m+1) 

11 


[In 

\in-h: 


rrm tt 

^22 ■ ^11 


(m + 1)] 


H 


21 


rr(m+l) 
22 ’ 11 


I -1 


H 


22 


Tjr(m+1) ^ 
■ 12 5 


tr(m) tt 

■^21 -^21 


(m+1) 


T rrm rj(m+l) 

J- ^22 ’ ^ 11 




(m+1) 




(m+1) 

21 


I -1 


'■ 22 


H 


215 

r(m+l) 
• 11 


-1 


rrm rr(m + l) 
’ 22 ' 12 


In the same way the Stiffness matrix relating the field and the linear form in the 
positions Zm-i and Zr can be described as the Stiffness matrix of the structure formed 
by the layers m, m + 1 , ..., /x: 


(4.12) 


A(m : Zm-i) 
A{fi : Zr) 


= 


F(m : Zm-i) 
F(^ : Zr) 


and its composition rule in terms of and E"^ is given by: 


(4.13) 




I7i(m) _ iTim 

£/ll — 


-^12 “ 


jp{m) _ p,(m+l) 

-^21 “ -^21 


j^{m) _ j^(m+l) 
■^22 ~ ^22 


E 


(m+1) 


11 


E 


22J -^215 


-C/9 


-,(m + l) 


-i(m+l) 


^(m+1) 


>22 


r- 






^22 


-1 


^(m-n) 


Analogously, the composition rule for the Scattering matrix S'(R;m), can be 
expressed in terms of the matrices S'(R;m+ 1) and ^(m + l;m), each one defined 
in agreement with (4.3) for the interfaces placed between the domains m and R, 
between m+1 and R and between m and m+1, respectively. This composition rule 
can be expressed by means of the product denoted by in the form: 


(4.14) iS'(R; m) = S'(R; m + 1) S'(m + 1; m) . 

Given three matrices Z, Y and X of order 2N subdivided in their NxN partitions 
the product @ expressing Z = Y @ X is defined in [4] by means of the composition 
rule: 

Zn = Xn + Xi 2 • Yn • [ 1 ^ - X22 ■ Yn]-' • X21; 

Z 12 = X 12 • Yi2 + X 12 • Yn • [In - X 22 • Yn]”^ • X 22 • Y^; 

Z 21 = Y 2 i-[Jw-X 22 -Yn]“'-X 2 i; 

(4.15) Z 22 = Y 22 + Y 21 • [In - X 22 • Yn]'" • X 22 • Y 12 . 

We must note that the composition rules (4.11) and (4.15) include the inverses 
I^Jat — H ^2 ' and [In — S' 22 (m + 1; m) ■ S'ii(R; m + 1)]”^ respectively, which 














TRANSFER MATRIX NUMERICAL STABILITY AND THE Cd PROBLEM 


13 


are regular even when the thickness of the layer or of the multilayer goes to infinity or 

-I -1 


(m+l) 


-C'9 


which is 


to zero. The composition rule (4.13) includes the term E 
regular when the thickness of the layer or of the multilayer goes to infinity. For very 
small thicknesses this composition rule will lead to the accumulation of the roundoff 
errors. 


5. General formulation of some typical boundary problems. Numerical 
examples. The boundary problems can be formulated in terms of the T, hybrid, scat¬ 
tering or stiffness matrices. We consider a system formed by three domains L —M —R. 
The internal domain M can be formed by one or several homogeneous layers in whose 
case the matrix of the structure M must be obtained through composition rules (see 
Section 4.2). In the external domains L and R we can have different media and even 
the vacuum. Depending on the problem under study we shall employ different bound¬ 
ary conditions at the interface L|M with coordinate zi and at M|R with coordinate Zr- 
We denote by F{i/r), A(£/r) the field and the associated linear form at the coordinate 
zijZr- In all the cases here considered we avoid to use submatrices which can exhibit 
numerical instabilities when d ^ oo, as it happens for oi' ^12 ■ 

5.1. Escape problem. We shall study the escape problem in a system formed 
by three media L — M — R having full matching conditions (FMC) at the interface 
L|M with coordinate Z(, and at the interface M|R with coordinate Zr- In the scape 
problem we shall have only outgoing waves in M. Applying the continuity conditions 
at the interface we can write: 


(5.1) 


F(f)- 
A(r) + 


HirJ) 


A{e)- 
F(r) + 


The superindex ± denote the vectors related with the wave travelling in R/L 
towards the right/left. From the two matrix equations coming from (5.1) we can 
write: 


(5.2) 


f -In HnirJ) 0^ 

Ojv H2i{r,£) H22{rJ) -In 


m- 

A{e)- 
F(r) + 
A(r) + 


We can express the vectors appearing in the right-hand side of (5.2) in the form: 


(5.3) 


(5.4) 






ai(L)~ 

F(f)- 


^ Fii£)- . 

.. Fn{£)- 


A{£)- 


_ Ai(^)- . 

.. An(£)- 


aAf(L)~ 









ai(R) + 

F{r) + 


Fi(r)+ . 

. FN{r) + 


A(r) + 


Ai(r)+ . 

■ Ajv(r)+ 






aAr(R)'^ 


= LI(£)- ■ a(L)' 


= LI(r)+ ■ a(R)- 


where Fj{£) are LI solutions belonging to the L domain, evaluated at zg. and Fj(r) + 
are LI solutions belonging to the R domain, evaluated at Zr- The A-vector a(R)+ is 
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formed by the coefficients aj(R)+ from those waves travelling to the right at R and 
a(L)“ by the coefficients aj(L)~ from those waves travelling to the left at L. 

Then by using (5.3) and (5.4) we transform (5.2) into the secular system: 


(5.5) 


( Msn 

Msi2 \ ^ 

a(L)” 

V Ms21 

MS 22 j 

a(R)+ 


(5.6) 

Msii 

(5.7) 

Msi2 

(5.8) 

Ms21 

(5.9) 

Ms22 


[-In Hn{rJ) ]-Ll{ey 
[Hi2ir,i) 0 ^].LI(r)+, 
[ 0 ^ // 2 i(r,^) ] •LI(£)-, 

[ ff22(r,£) -/iv]-LI(r)' 


The problem eigenvalues are obtained from the secular equation Det[Ms] = 0. 
In terms of the Stiffness matrix we have: 


(5.10) 

Msii 

= [ -Eii(r,£) 

-In ] • LI(^)- 

(5.11) 

Msi2 

= [ -Ei2(r,£) 

0^ ] •LI(r)+, 

(5.12) 

Ms21 

= [ E2i{r,£) 

0^ ] •LI(£)-, 

(5.13) 

Ms22 

= [ E22{r,£) 

-In ] • LI(r) + 


As a numerical example we use the secular equation in terms of the hybrid matrix 
H{r, £) to obtain the velocities of shear horizontal (SH) acoustic waves in fx piezoelec¬ 
tric multilayers systems. These curves were obtained in Ref. [3] by using the singular 
value decomposition (SVD) method together with a variant of the Global Matrix 
Method (GMM) as an alternative technique to avoid the numerical instabilities found 
by the authors. 

The piezoelectric systems studied there, are formed by two different materials, A 
(PZT4) and B (PZT5A), and have different layer configurations: n = 3 (ABA), n = 5 
(ABABA), n = 7 (ABABABA) and n = 9 (ABABABABA). All these systems have 
the L — M — R structure, with L = R = A. The external domains are semi-infinite 
and to obtain confined modes it was assumed that there are no ingoing waves in the 
inner region M, whereas the outgoing waves are evanescent. It is then clear that 
this problem can be studied as a particular case of the scape problem considered in 
this section. In order to get evanescent waves the eigenvalues kj appearing in the 
exponential terms of these waves were assumed to be pure imaginary. 

Except for n = 3, the hybrid matrix H (r, £) in the inner region M was obtained 
by means of the composition rule (4.11). To solve this problem it was necessary 
to transform the original system of two equations of motion [3] in a matrix Sturm- 
Liouville system (2.1) with N=2. In this problem F(z) has two components, the 
transverse displacement u and the electric potential (/). The 2 axis is oriented in the 
direction normal to the multilayer interfaces in such a way that it coincides with the 
y axis in the scheme of Figure 1 in Ref. [3]. The x axis coincides in both cases. 

The quadratic eigenvalues problem solution (QEP, Section 2.1) for one layer is: 
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(5.14) 

(5.15) 

(5.16) 

(5.17) 


ki 


= -iKx 



k2 = -ki 


k?. = 


\ 


— Ki + 


(C44 + ^) 
£11 


^4 = -^3 


where Vg is the velocity of the surface wave we are studying, whereas v = \ (C44 H - )/p 

Y £11 

is the SH wave velocity. The material parameters of the layer needed for this study 
are the mass density p, the elastic constant C 44 , the piezoelectric constant 615 and the 
dielectric constant en. The eigenfunctions can be chosen in the form: 


(5.18) Fj{z) = Fjo ^ ( 1 ) e*''"J = 1,2 

and: 

(5.19) Fj{z) = F,o ^ ^ = 3 ^ 4 ^ 

fca and fe 4 must be pure imaginary to obtain evanescent outgoing waves. As these 
waves travel in material A (PZT4) layers the expression (5.16) shows that this happens 
for Vg < VA- It is also possible to obtain confined modes when there are layers in the 
domain M with k^ and ^4 real. This is only possible in material B (PZT5A) layers 
when VB < Vg. 

The hybrid matrix of an independent layer was obtained by a method analogous 
to that employed in [4] to get the expression (2.8). For the H matrix we have: 


(5.20) 

H{z,zo)=V^^{z 

^o) • [U^^(^ 

^ 0 )] ^ 

(5.21) 

ttFA/ \ Fi(zo) 

U 

F 2 ( 2 o) 

A 2 (z) ... 

F2Ar(2:o) 

A2Ar(z) 

(5.22) 

^ \ Ai( 2 :o) 

A 2 ( 2 :o) ... 

F2{z) . . . 

A27v(^:o) 

F2Ar(^) 


The secular matrix Ms was obtained from the expressions (5. 6 -5.9) and then we 
obtained the values of the surface wave velocities zeroing the secular determinant at 
different frequency values. Table 1 shows the values obtained in our calculation, those 
obtained in [3] together with the corresponding frequencies. The values in [3] were 
obtained by using the (SVD) method and a Global Matrix of order 4:{Nb — 1) x 4(Ab;, — 
1), Nb being the number of layers in the structure. Thus for Nb = 3 the matrix would 
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No. of 

UJ 

MG and SVD 

Hs 

layers 

(MHz) 

vs (m/s) 

vs (m/s) 

3 

123.1 

2324 

2324.08 


357.1 

2286 

2285.94 

5 

123.1 

2313.6/ 2340.5 

2313.9/ 2340.8 


279.1 

2292.6/ 2294.5 

2292.5/ 2294.7 


318.1 

2343.1/ 2350.7 

2344.1/ 2351 


396.1 

2330.3/ 2333.6 

2330.3/ 2333.8 

9 

20 

2339 

2339.4 


80 

2314/ 2335 

2314/ 2335.2 


Table 1 


Comparison between the surface wave velocity values for different frequency values obtained by 
two different theoretical methods: (GM) Global Matrix Method and (SVD) Singular Value Decom¬ 
position method. (HsJ Hybrid compliance-stiffness Matrix Method. 



Fig. 2. Surface wave velocity values for different frequency values of the n = 3 and n = 9 systems 


be (8 X 8 ), whereas for Np = 9 the matrix would be (32 x 32). The hybrid matrix 
employed in our calculations is of order (4x4). The good agreement of both sets 
of velocity values shows the capability of the hybrid matrix method to avoid the VLd 
problem with lower computational and formal requirements when compared with the 
Global Matrix method. 

Figure 2 shows the values of the surface wave velocity for the corresponding 
frequency values for the three and nine layer systems coming from our calculations. 
We can observe two bands, the first of the even modes and the first of the odd modes, 
together with the convergence of the modes of the system A^=9 towards those of the 
system iV=3 when the frequency is increased. This behaviour is present in the curves 
given in [3]. 
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5.2. Periodic systems. Let us consider a periodic system along the 2 direction 
with arbitrary period d. This could be a periodic bulk crystal or a superlattice. The 
matrices H{z + d,z), E{z + d,z) transfer along a given period. The Bloch-Floquet 
conditions are satisfied for both F( 2 :) and A( 2 :), in such a way that F( 2 :+(i) = F{z)-e^‘^‘^ 
and A{z + d) = A{z) ■ 

We can write this in terms of the hybrid matrix as: 


(5.23) 


F(z) 
A{z + d) 


H{z + d,z) 


A{z) 
F(z + d) 


The Bloch-Floquet conditions for F(z) and A( 2 ;) in (5.23) lead to: 


(5.24) A(z) = ■ [I - ■ F(z) 

(5.25) A(z) = [I - • H22 ■ F{zy,. 

We write H instead oi H{z + d, z) to simplify. The secular system is obtained 
from expressions equations (5.24) and (5.25): 

(5.26) {[7 - 772ie-*«‘'] • H 22 - ■ [I - 77i2e*«^]} • F(z) = 0^ . 

It will have nontrivial solutions if: 

(5.27) Det { [7 - 772ie-*«''] ■ H 22 - ■ [7 - 77i2e*«'^]} = 0. 

This equation gives a dispersion relation in terms of the 77 matrix elements for 
any N. 

Following the same procedure with E{z + d,z) we obtain the following secular 
system: 


(5.28) {[Ell + Ei 2 e^‘i‘^] - [E 2 ie-^‘‘‘^ + E 22 ]} • F(2) = Ojv , 

and the dispersion relation: 


(5.29) Det {[Ell + - [7;2ie“*«‘^ + E 22 ]} = 0 

We assume that in our periodic system the inner domain M (containing one or 
several homogeneous layers) coincides with the period d. Now we shall pose the 
problem in terms of the Scattering matrix S'(R; L). 

For the external domains L and R we have: 


(5.30) 


F(L : Zi) 
A(L : Zi) 


Q{L : ze) ■ 


a+(L) 
a (L) 


(5.31) 


F(R : Zr) 
A(R : Zr) 


(5(R : Zr) ■ 


a+(R) 

a-(R) 
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From now on we shall employ QL instead of Q(L : zg) and QR instead of Q(R : Zr) 
to simplify the notation. Usually a reduced base in zg is employed to obtain QL and 
a reduced base in Zr is used to obtain QR. From these matrices we can obtain the 
matrix /^(R, L) by means of (4.5) and then from (4.4) we can obtain S'(R; L). 

From the Bloch-Floquet condition we obtain: 


(5.32) 


F(R : Zr) 


F(L : 

A(R : Zr) 


A(L : zg) 


Combining (5.30), (5.31) and (5.32) with the expression (4.3) defining the Scat¬ 
tering matrix we can write the following expressions: 


a+ (L) = [QRii . S21 - QLii - QL12 • 5ii • 

(5.33) [QL12 • S'12 - QRii • S22 - QR12] ■ a"(R). 

a+(L) = [QR21 • ^21 - QL21 - QL22 • Sn ■ 

(5.34) [QL22 ■ 5'12 — QR21 • S22 ^ QI422] ■ a (R). 

Subtracting these equations we arrive to the secular system: 


Otv = Ms ■ a (R), 

and from it we obtain the secular determinant: 

Det {[QRii ■ S21 - QLn - QL12 • 5 ii • 

[QL12 • 5 i 2 - QRii ■ S22 - QR12] 

- [QR21 ■ 521 - QL21 - QL22 ■ 5 ii • 

(5.35) [QL22 • 5 i 2 — QR21 • S22 ^ QR22]} = 0 . 

We note that the equations (5.27), (5.29) and (5.35) are given in terms of matrix 
blocks that can overcome the numerical instability known as fid problem. This is not 
the case for the secular equation in terms of T( 2 : -I- d, z): 


(5.36) Det[T(z + d,z)-I = 0 , 

We shall consider now as an example the motion of electrons in a periodic one¬ 
dimensional potential such as that of a superlattice formed by barriers of B material 
with effective mass ms, thickness b and height Vq and wells of A material with effective 
mass ruA and thickness a. In this case the equations (5.27), (5.29), (5.35) and (5.36) 
are given by: 
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(5.37) 2 cos{qd)Hi2 = 1 - H 11 H 22 + HI 2 

(5.38) 2 cos((7(i)£^i2 = E 22 ^ En 

(5.39) 2cos{qd)Si2 = {S 21 S 12 - SnS 22 ) + 1 

(5.40) cos(gd) = i(rn+T 22 ), 
where 

(5.41) kA = 

(5.42) 

(5.43) kB = \l^{E-Vo). 

We used the sin k(z — ztjZy) and cos k{z — znjZy) base in the L/R domain to obtain 
(5.39). When we use the matrix elements of T for this problem in the period d = a + 5 
in (5.40) we arrive to the well known Kronig-Penney equation [31]. The expressions 
(5.37)-(5.39) are variations of this equation if we notice that the matrices iT, E and 
S can be calculated from their relations with T. 

Expressions (5.37)-(5.39) are variations of (5.40) to calculate the system energy 
levels for any barrier width b. When the barrier thickness b ^ co (limit of isolated 
symmetric rectangular wells) the secular equation in terms of T diverges. On the 
other hand its variations lead directly to the well known transcendental equations 
giving the energy levels for even and odd states of a symmetric rectangular well of 
width a and depth Vq. 

After some algebra it can be shown that equation (5.37) coincides with the equa¬ 
tion (32) in [32] for the Kronig-Penney equation. In the same way it coincides with the 
equation (20) of Ref. [33]. Refs. [32, 33] give results for the Kronig-Penney equation 
to avoid the rid problem. 

6. Conclusions. In the general framework of N equation systems of the Sturm- 
Liouville matrix kind with piecewise constant coefficients we have shown that there are 
transfer matrix variants with dimensions independent of the number of layers in the 
structure which can avoid the numerical instabilities present in the ATM. The hybrid 
compliance-stiffness matrix and the scattering matrix can avoid the so called rid 
problem, being numerically stable independently of how big or small be the thicknesses 
in the multilayer structure. The Stiffness matrix and its inverse the compliance matrix 
are numerically stable for big thicknesses of the layers or of the multilayer structure. 
On the other hand, in the case of very small layer thicknesses these two matrices can 
exhibit the nd problem due to the roundoff errors accumulation.For zero thicknesses 
both matrices exhibit a numerical singularity (overflow). 

Given the big variety of boundary problems which can be studied with these 
numerically stable variants of the ATM and the generality and ubiquity of the matrix 
Sturm-Liouville system, the results obtained here can be applied to the study of 
various elementary excitations in multilayer systems. 

The relations between the different matrices studied here has proven to be an 
useful instrument in the study of the numerical stability of transfer matrices. With 
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this technique it was possible to show analytically the capability of some of these 
variants of the transfer matrix to avoid the numerical degradation leading to the fid 
problem. 

In recent years some methods able to deal with systems having inhomogeneous 
coefficients have been developed. We present in Appendix C the link of the N equation 
systems of the Sturm-Liouville matrix kind to the corresponding differential forms of 
those problems. 

Appendix A. Example of a matrix subset with a similar behaviour from 
the numerical point of view. 

Let us denote by A, y, Z and R four matrices in whose definition enter the 
vectors F( 2 :), F(zo), A{z) and A( 2 :o), as for example: 


A(zo) 

F(z) 

= X 

F{zo) 

A{z) 


F(z) 

A(zo) 

= Y 

F{zo) 

A{z) 

A(zo) 

F(^) 

= z 

A{z) 

F(^o) 


F(z) 

A( 2 :o) 

= R 

A{z) 

F(^o) 


If we take as the reference matrix any one of them it can be shown that one of 
the remaining matrices is obtained by permutations among them of the vectors in the 
right-hand side of the reference matrix. A second one is obtained by following this 
method among the vectors in the left-hand side of the reference matrix. Finally the 
third one is obtained with both permutations. The relations between these matrices 
can be resumed as: 


(A.2) (A)n = (y)2i = iR)22 = iZ)i2 

(A.3) (A)i2 = (y)22 = {R)21 = {Z)ll 

(A.4) (A)21 = (y)ll = (i^)l2 = {Z)22 

(A.5) (A)22 = (y)l2 = {R)ll = {Z)21 

These relations show that the matrices X, Y, Z and R will have a similar 
behaviour from the numerical point of view. 

Appeudix B. Matrix T partitions of order N. 

Starting with the expression: 


(B.l) T{d) 


^Oat ^02iV 


nfc;^(d) 0 


Fqa^ Fo2Af 

■^Oat ■^02A/^ 


0 nfe2w(d) 


Aqat ■^02A/^ 


we have: 


Til = Foat ■ Flfej, (d) • + F02JV ■ Flfejiv ('^)7i2 ■ 

Ti2 = Foat ■ FIfej, (d) • 721 + Fqjat ■ FIfc2jv(d)722 • 

T 21 = Xqj^ • (d) • 7 ii + Ao 2 Ar ■ FIfc 2 Ar ('^) 7 i 2 • 

T 22 = Aojv • nfej^(d) • 7 ^ 1 ^ + Ao 2 „ ■ Ilk 2 N{d)j 22 ■ 


(B.2) 
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Til ~ [^^OiV ^^02iv ' ■^02N ' -^Ow]- 
■712 ~ [^^02iv ■ -^Oat ’ ■^02Ar]- 

T 2 I ~ [-^On -^02N ■ ^02Ar ' ^On]- 

(B- 3 ) 722 = [Ao 2„ - Ao„ • ■ Fo2a,]- 


Appendix C. Sturm-Liouville matrix form for inhomogeneous media. 

The matrix Sturm-Liouville equation 


(C.l) 
_d 
dz 


B{z) 


dF{z) 

dz 


+ P{z)-F{z) 


+ Y{z) 


dF{z) 

dz 


+ W{z)-F{z) = U^)-F{z) = Qn>.i 


can be written as: 


(C.2) 

^ F{z) 
dz a{z) 


-B{z)-'^P{z) 

B{z)-^ \ 

F{z) 

Y{z) ■ Biz)-Y p{z) - W{z) 

-Y{z)-B{z)-^ ) ■ 

A{z) 


Here A{z) 


B{z) 


dF{z) 

dz 


+ P{z) ■ F{z) is the SLM operator matrix differential 


form. 

The equation (C.2) is the link with the first order differential equations systems 
given in eq.(2.10) of [16], eq.(3.3) of [17], eq.(2) and (A.6) of [18], eq.(3) of [19] and 
eq.(8) of [20]. 

These equations cover different inhomogeneous systems. 


C.l. Radially inhomogeneous cylindrically anisotropic systems. This is 
the case considered in [16]. In this work the mass density and the elements of the 
stiffness tensor depend only on the radial coordinate r. It is then possible to write: 


(C.3) 




to obtain: 


d 

dr 


^diCU) 

rQ —- + (Rk + zK^rP) 

dr 


kRr + iKzrP^ 


djCU) 

dr 


(C.4) 


+ - 


fcTfc + iKzr{kS + S^k) + (m^r)^(M - IpLoyni) (CU) = 0. 


Here C is a normalization constant. Equation (C.4) is of kind (C.l) for a cylin¬ 
drical elastic material radially inhomogeneous. Here r plays the role of z in the planar 
systems. With the properties imposed on fc, Q, P, R, S, T and M in [16] we have 
Y = -P\ B = B^ &ndW = Wl 

In this case the linear differential form in (C.4) is A = rtr (where t^) is the stress 
radial component and from eq.(2.7) from [16] we obtain A(r) = C'rY^"^(r). We see 
also that F{r) = CU^^\r). By identifying the B, P, Y and W matrices in (C.4) 
and substitution in (C.2) we arrive to eq.(2.10) of [16], eq.(3.3) of [17] and eq.(A.6) 
of [18]: 
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(C.5) 


= -G(r)T7(r)(”); 
dr r 


(C.6) 


= C 


U^^\r) 

irT^^Hr) 


(C.7) -G(r) 

r 


Y{r) ■ -B(r)-i • P(r) - W(r) 


-y(r) . ; • 


C.2. Shear-horizontal elastic waves in phononic crystals formed by in- 
homogenoeus anisotropic materials. Cartesian coordinates. This is studied 
in [20] where the displacement u depends on xi and X 2 ^ but after expanding xi in 
plane waves they obtain the following ordinary differential equation in X 2 


(C.8) — {di + iKi)(p(5i + iKi)u) -I- d2{fJ-d2u) = —pto^u. 

We can then identify F{x 2 ) = u{x 2 ) and ^( 0 : 2 ) = pd 2 U, P = 0,Y = 0, B = p~^ 
and W = puP — (^1 + iKi)p{di +iKi). After substitution of these expressions in (C.2) 
W acts on the displacement u{x 2 ) and we obtain eq.(8) of [20] which is essentially 
the same than eq.(3) of [19]. 

We have seen that in all these cases involving inhomogeneous elastic anisotropic 
media we can put the matrix Sturm-Liouville in the (C.2) form. Then it would be 
possible to apply the stable integration methods of [16, 17, 18] for cylindrical geometry 
and those of [19, 20] for layered systems. 

Acknowledgments. We thank the Associate Editor and the referees for valuable 
comments and useful suggestions. 
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